Using GRASS GIS in Jupyter Notebooks:
An Introduction to grass.jupyter

Caitlin Haedrich, Vaclav Petras, Anna Petrasova

grass_jupyter

NCSU GeoForAll Lab at the Center for Geospatial Analytics
North Carolina State University

Caitlin Haedrich

  • PhD student GeoForAll Lab at North Carolina State University
  • Advised by Helena Mitasova (and Anna Petrasova and Vashek Petras)
  • Interested in geovisualization and increasing accessibilty of GIS tools
  • Working with integrating GRASS GIS with Jupyter through Google Summer of Code 2021 and Mini Project Grant from GRASS GIS
  • TA for graduate-level Geospatial Computation and Simulation

geoforall

ncstate

GRASS GIS

  • Open-source GIS software with over 500 tools and 400 addons
  • Interfaces: graphical user interface, command line, Python, C
  • 3rd Party Interfaces: R, QGIS, actinia (REST API)

GUI

grassgis

Project Jupyter

“The Jupyter Notebook is an open-source web application that allows you to create and share documents that contain live code, equations, visualizations and narrative text.” – jupyter.org

Jupyter

In [1]:
import sys

v = sys.version_info
print(f"We are using Python {v.major}.{v.minor}.{v.micro}")

We are using Python 3.9.5

The output can also be interactive.

In [2]:
from ipywidgets import interact

def f(x):
    return x

interact(f, x=10);

interactive(children=(IntSlider(value=10, description='x', max=30, min=-10), Output()), _dom_classes=('widget-…

Teaching Geospatial Analytics with Jupyter

GIS714: geospatial computation and simulation: graduate level course required for Geospatial Analytics PhD students

oldGIS714

grass.jupyter

  • Released with main GRASS GIS distribution starting in version 8.2
  • Python package that makes GRASS more accessible in Jupyter by providing:
    • Data management functions
    • Visualizations classes for viewing and interacting with GRASS through in-line figures
  • Created with syntax consistent with GRASS's existing Python API and command line interface


You can try it out in Binder! (also linked from the GRASS GIS GitHub README) Binder

Getting Started

In [3]:
# Import Python standard library and IPython packages we need.
import subprocess
import sys

# Ask GRASS GIS where its Python packages are.
sys.path.append(
    subprocess.check_output(["grass80", "--config", "python_path"], text=True, shell=True).strip()
)
In [4]:
# Import the GRASS GIS packages we need.
import grass.script as gs
import grass.jupyter as gj

# Start GRASS Session
session = gj.init("data", "nc_spm_08_grass7", "PERMANENT")
In [17]:
# Set region to the high resolution study area
gs.run_command("g.region", region="rural_1m")

And then we are ready to begin using GRASS!

We can use the command line interface for GRASS by using the cell magic ! (for individual lines) or %%bash (applies to entire cell) which sends the code to the terminal.

In [5]:
!g.list type=raster -m -t

raster/aspect@PERMANENT
raster/basin_50K@PERMANENT
raster/boundary_county_500m@PERMANENT
raster/cfactorbare_1m@PERMANENT
raster/cfactorgrow_1m@PERMANENT
raster/el_D782_6m@PERMANENT
raster/el_D783_6m@PERMANENT
raster/el_D792_6m@PERMANENT
raster/el_D793_6m@PERMANENT
raster/elev_lid792_1m@PERMANENT
raster/elev_ned_30m@PERMANENT
raster/elev_srtm_30m@PERMANENT
raster/elev_state_500m@PERMANENT
raster/elevation@PERMANENT
raster/elevation_shade@PERMANENT
raster/facility@PERMANENT
raster/geology_30m@PERMANENT
raster/lakes@PERMANENT
raster/landclass96@PERMANENT
raster/landcover_1m@PERMANENT
raster/landuse96_28m@PERMANENT
raster/lsat7_2002_10@PERMANENT
raster/lsat7_2002_20@PERMANENT
raster/lsat7_2002_30@PERMANENT
raster/lsat7_2002_40@PERMANENT
raster/lsat7_2002_50@PERMANENT
raster/lsat7_2002_61@PERMANENT
raster/lsat7_2002_62@PERMANENT
raster/lsat7_2002_70@PERMANENT
raster/lsat7_2002_80@PERMANENT
raster/ncmask_500m@PERMANENT
raster/ortho_2001_t792_1m@PERMANENT
raster/roadsmajor@PERMANENT
raster/slope@PERMANENT
raster/soilsID@PERMANENT
raster/soils_Kfactor@PERMANENT
raster/streams_derived@PERMANENT
raster/towns@PERMANENT
raster/urban@PERMANENT
raster/zipcodes@PERMANENT
raster/zipcodes_dbl@PERMANENT

Or we can use the Python API for GRASS GIS. Since grass.jupyter is written in Python, we'll use the GRASS Python API from here on out.

In [6]:
gs.read_command("g.list", type="vector", flags="mt")
Out[6]:
'vector/P079214@PERMANENT\r\nvector/P079215@PERMANENT\r\nvector/P079218@PERMANENT\r\nvector/P079219@PERMANENT\r\nvector/boundary_county@PERMANENT\r\nvector/boundary_municp@PERMANENT\r\nvector/bridges@PERMANENT\r\nvector/busroute1@PERMANENT\r\nvector/busroute11@PERMANENT\r\nvector/busroute6@PERMANENT\r\nvector/busroute_a@PERMANENT\r\nvector/busroutesall@PERMANENT\r\nvector/busstopsall@PERMANENT\r\nvector/census_wake2000@PERMANENT\r\nvector/censusblk_swwake@PERMANENT\r\nvector/comm_colleges@PERMANENT\r\nvector/elev_lid792_bepts@PERMANENT\r\nvector/elev_lid792_cont1m@PERMANENT\r\nvector/elev_lid792_randpts@PERMANENT\r\nvector/elev_lidrural_mrpts@PERMANENT\r\nvector/elev_lidrural_mrptsft@PERMANENT\r\nvector/elev_ned10m_cont10m@PERMANENT\r\nvector/firestations@PERMANENT\r\nvector/geodetic_pts@PERMANENT\r\nvector/geodetic_swwake_pts@PERMANENT\r\nvector/geology@PERMANENT\r\nvector/geonames_NC@PERMANENT\r\nvector/geonames_wake@PERMANENT\r\nvector/hospitals@PERMANENT\r\nvector/lakes@PERMANENT\r\nvector/nc_state@PERMANENT\r\nvector/overpasses@PERMANENT\r\nvector/poi_names_wake@PERMANENT\r\nvector/precip_30ynormals@PERMANENT\r\nvector/precip_30ynormals_3d@PERMANENT\r\nvector/railroads@PERMANENT\r\nvector/roadsmajor@PERMANENT\r\nvector/schools_wake@PERMANENT\r\nvector/soils_general@PERMANENT\r\nvector/soils_wake@PERMANENT\r\nvector/streams@PERMANENT\r\nvector/streets_wake@PERMANENT\r\nvector/swwake_10m@PERMANENT\r\nvector/urbanarea@PERMANENT\r\nvector/usgsgages@PERMANENT\r\nvector/zipcodes_wake@PERMANENT\r\n'

Visualizing GRASS data

In [ ]:
In [ ]:
In [ ]:

grass.jupyter Syntax

  1. Create Instance

  2. Add layers:

    d.rast -> map.d_rast(), map3D.d_rast

  1. show and save
In [ ]:

folium Integration

folium...

  • is a popular Python library for creating Leaflet.js maps
  • creates interactive HTML maps that can be displayed in-line
  • has built-in tilesets from OpenStreetMap, Mapbox, and Stamen, and supports custom tilesets

foliumLogo

Moving data from GRASS GIS location projection to folium is a challenge

folium is projected in Web Mercator (EPSG:3857)

However, any coordinates (i.e. vector data) need to be specified in degrees of latitude and longitude (WGS84, EPSG:4326)

folium reprojects latitude and longitude to Web Mercator internally



We pass data to folium by reprojecting to temporary locations

grass_folium

In [7]:
import folium 

# Create figure
fig = folium.Figure(width=600, height=400)

# Create a map to add to the figure later
m = folium.Map(tiles="Stamen Terrain", location=[35.761168,-78.668271], zoom_start=13)

# Create and add elevation layer to map
gj.Raster("elevation", opacity=0.5).add_to(m)

# Do some cool folium stuff!
# Like make a tooltip
tooltip = "Click me!"
# and add a marker
folium.Marker(
    [35.781608,-78.675800], popup="<i>Center For Geospatial Analytics</i>", tooltip=tooltip
).add_to(m)

# and a circle
folium.Circle(
    radius=120,
    location=[35.769781,-78.663160],
    popup="Great Picnic Area",
    color="crimson",
    fill=False,
    tooltip=tooltip
).add_to(m)

# Add the map to the figure
fig.add_child(m)

# Display figure
fig
Out[7]:

The InteractiveMap class allows users unfamiliar with folium to produce maps easily.

In [8]:
# Create Interactive Map
fig = gj.InteractiveMap(width = 600)
# Add raster, vector and layer control to map
fig.add_raster("elevation", title="Elevation Raster", opacity=0.8)
#fig.add_vector("roadsmajor")
fig.add_layer_control(position = "bottomright")
# Display map
fig.show()
Out[8]:

Teaching Geospatial Computation and Simulation with grass.jupyter

GIS714: Geospatial Computation and Simulation, spring 2022,


Course Repo on GitHub

Binder

In [9]:
# Make a new mapset for this assignment
gs.run_command("g.mapset", mapset="HW3_water_simulation", location="nc_spm_08_grass7", flags="c")

Create a map of flooded area¶

We create a map of flooded area with r.lake (GRASS manual for r.lake) by providing a water level and a seed point:

In [21]:
gs.run_command("r.lake", elevation="elev_lid792_1m", water_level=113.5, lake="flood1", coordinates="638728,220278")

# See results
flood1 = gj.Map(use_region=True)
flood1.d_rast(map="elev_lid792_1m")
flood1.d_rast(map="flood1")
# Display map
flood1.show()
Out[21]:
Question 4¶

Increase water level to 113.7m and 114.0m and create flooded area maps at these two levels.

In [22]:
#### Your Answer Here

Estimating inundation extent using HAND methodology¶

We will use two GRASS addons, r.stream.distance and r.lake.series, to estimate inundation with Height Above Nearest Drainage methodology (A.D. Nobre, 2011). We need to install them first:

In [ ]:
#!g.extension r.stream.distance
#!g.extension r.lake.series

For this section, we will change our computation region to elevation which is a larger study area than we used above. Because we are using a new region (and also a higher threshold of 100,000), we need to run r.watershed again to compute the flow accumulation, drainage and streams. We convert the streams to vector for better visualization.

In [24]:
gs.run_command("g.region", raster="elevation")
gs.run_command("r.watershed", elevation="elevation", accumulation="flowacc", drainage="drainage", stream="streams_100k", threshold=100000)
gs.run_command("r.to.vect", input="streams_100k", output="streams_100k", type="line")

Now we use r.stream.distance with output parameter difference to compute new raster where each cell is the elevation difference between the cell and the the cell on the stream where the cell drains. This is our HAND terrain model.

In [ ]:
gs.run_command("r.stream.distance", stream_rast="streams_100k", direction="drainage", elevation="elevation", method="downstream", difference="above_stream")

With r.lake.series, we can create a series of inundation maps with rising water levels. r.lake.series creates a space-time dataset. We can use temporal modules to further work with the data...

In [ ]:
gs.run_command("r.lake.series", elevation="above_stream", start_water_level=0, end_water_level=5, water_level_step=0.5, 
               output="inundation", seed_raster="streams")

... or visualize the flood with TimeSeriesMap.

In [ ]:
flood_map = gj.TimeSeriesMap(use_region=True)
flood_map.add_raster_series("inundation")
flood_map.d_legend(color="black", at=(10,40,2,6)) #Add legend
flood_map.d_barscale()
flood.show()

Course Reflections

Notebook format was generally well received:

  • Several students noted they enjoyed being able to experiment with GRASS tools by varying parameters and visually observing the effect while working through the tutorial sections of the notebooks.
  • By hosting the notebooks on GitHub and including a link to Binder, students could easily access and run the notebooks without installing and configuring software on their personal computers.
    • This also gave the students exposure to working with Git and GitHub, an important tool for open science and collaborative research.

There were several areas of the notebook deployment that needed improvement or that posed a difficulty to the class:

  • First, running the notebooks locally was difficult for many students. Due to the required configuration of path variables on Windows operating system and installation of required packages, many students switched to Binder right away.
  • For students that did get the notebooks running locally, we often encountered issues that were specific to their version of GRASS GIS or operating system.
  • We also found that Binder was not a reliable host; occasionally, there were too many users accessing Binder’s computational nodes and the notebooks failed. In the future, we may use a JupyterHub hosted by NC State University to circumvent such technical issues.
  • Secondly, the notebooks assume users have basic Python abilities, are familiar with the fundamentals of geospatial analysis (such as common data formats and projections) and have some experience using Jupyter Notebooks. For students without any of these prerequisite skills, the course could be challenging even with the removal of the command line interface and GUI workflows present in the previous tutorial format.

More Resources

grass.jupyter Manual Page

FOSS4G '22 GRASS GIS Workshop Materials (Anna Petrasova)

GIS714: Advanced Geospatial Computation and Simulation Course Repository

Example Notebooks in GRASS GIS's GitHub Repository (can be run with binder through link provided in main repo README)

QRcode
Find this presentation on GitHub

Future Work¶

Thank you to...¶

Helena Mitasova

Stefan Blumentrath

Vero Andreo

GIS714 Class Spring 2021

… and the GRASS community, NC State Center for Geospatial Analytics, Google Summer of Code